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Abstract 

A single mammalian cell includes an order of 10 4 —10 5 mRNA molecules and as many as 10 5 —10 6 ribosomes. 
Large-scale simultaneous mRNA translation and the resulting competition for the available ribosomes has important 
implications to the cell’s functioning and evolution. Developing a better understanding of the intricate correlations 
between these simultaneous processes, rather than focusing on the translation of a single isolated transcript, should 
help in gaining a better understanding of mRNA translation regulation and the way elongation rates affect organismal 
fitness. A model of simultaneous translation is specifically important when dealing with highly expressed genes, 
as these consume more resources. In addition, such a model can lead to more accurate predictions that are needed 
in the interconnection of translational modules in synthetic biology. 

We develop and analyze a general model for large-scale simultaneous mRNA translation and competition for 
ribosomes. This is based on combining several ribosome flow models (RFMs) interconnected via a pool of free 
ribosomes. We prove that the compound system always converges to a steady-state and that it always entrains 
or phase locks to periodically time-varying transition rates in any of the mRNA molecules. We use this model 
to explore the interactions between the various mRNA molecules and ribosomes at steady-state. We show that 
increasing the length of an mRNA molecule decreases the production rate of all the mRNAs. Increasing any of the 
codon translation rates in a specific mRNA molecule yields a local effect: an increase in the translation rate of this 
mRNA, and also a global effect: the translation rates in the other mRNA molecules all increase or all decrease. 
These results suggest that the effect of codon decoding rates of endogenous and heterologous mRNAs on protein 
production is more complicated than previously thought. 

Index Terms 

mRNA translation, competition for resources, systems biology, monotone dynamical systems, first integral, 
entrainment, synthetic biology, context-dependence in mRNA translation, heterologous gene expression. 


I. Introduction 

Various processes in the cell utilize the same finite pool of available resources. This means that the 
processes actually compete for these resources, leading to an indirect coupling between the processes. 
This is particularly relevant when many identical intracellular processes, all using the same resources, 
take place in parallel. 

Biological evidence suggests that the competition for RNA polymerase (RNAP) and ribosomes, and 
various transcription and translation factors, is a key factor in the cellular economy of gene expression. 
The limited availability of these resources is one of the reasons why the levels of genes, mRNA, and 
proteins produced in the cell do not necessarily correlate f[63ll . 021 . 1571 . If67l . If53l . [|69ll . Il29ll . 

It was estimated that in a yeast cell there are approximately 60, 000 mRNA molecules. These can be 
translated in parallel ll74l . ifTTl . with possibly many ribosomes scanning the same transcript concurrently. 
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Fig. 1. Illustration of simultaneous translation of mRNA chains (right) interconnected via a pool of free ribosomes (left). 


The number of ribosomes is limited (in a yeast cell it is approximately 240, 000) and this leads to a 
competition for ribosomes. For example, if more ribosomes bind to a certain mRNA molecule then the 
pool of free ribosomes in the cell is depleted, and this may lead to lower initiation rates in the other 
mRNAs (see Fig. []])• 

There is a growing interest in computational or mathematical models that take into account the com¬ 
petition for available resources in translation and/or transcription (see, for example, ll22ll . Il23l . If44l . ||9ll , 
1(131 . 10). One such model, that explicitly considers the movement of the ribosomes [RNAP] along the 
mRNA [DNA], is based on a set of asymmetric simple exclusion processes (ASEPs) interconnected to a 
pool of free ribosomes. ASEP is an important model from non-equilibrium statistical physics describing 
particles that hop randomly from one site to the next along an ordered lattice of sites, but only if the 
next site is empty. This form of “rough exclusion” models the fact that the particles cannot overtake one 
another. ASEP has been used to model and analyze numerous multiagent systems with local interactions 
including the flow of ribosomes along the mRNA molecule Il38l . ll58l . In this context, the lattice represents 
the mRNA molecule, and the particles are the ribosomes. For more on mathematical and computational 
models of translation, see the survey paper lf70ll . 

Ref. If24ll considered a closed system composed of a single ASEP connected to a pool (or reservoir) of 
“free” particles. The total number of particles is conserved. This is sometimes referred to as the parking 
garage problem , with the lattice modeling a road, the particles are cars, and the pool corresponds to a 
parking garage. Ref. lfT2l studied a similar system using domain wall theory. Ref. Iil3l (see also ll2Tfl ) 
considered a network composed of two ASEPs connected to a finite pool of particles. The analysis in 
these papers focuses on the phase diagram of the compound system with respect to certain parameters, 
and on how the phase of one ASEP affects the phase of the other ASEPs. These studies rely on the phase 
diagram of a single ASEP that is well-understood only in the case where all the transition rates inside the 
chain (the elongation rates) are equal. Thus, the network is typically composed of homogeneous ASEPs. 
Another model 0 combines non-homogenous ASEPs in order to study competition between multiple 
species of mRNA molecules for a pool of tRNA molecules. This study was based on the Saccharomyces 
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cerevisiae genome. However, in this case (and similar models, such as IfHlO analysis seems intractable 
and one must resort to simulations only. 

Our approach is based on the ribosome flow model (RFM) If52ll . This is a deterministic, continuous-time, 
synchronous model for translation that can be derived via the mean-field approximation of ASEP |f6]| . 
The RFM includes n state-variables describing the ribosomal density in n consecutive sites along the 
mRNA molecule, and n + 1 positive parameters: the initiation rate A 0 , and the elongation rate A* from 
site i to site i + 1, i = 1,..., n. The RFM has a unique equilibrium point e = e(A 0 ,..., A n ), and any 
trajectory emanating from a feasible initial condition converges to e ll42ll (see also If39l0 . This means that 
the system always converges to a steady-state ribosomal density that depends on the rates, but not on the 
initial condition. In particular, the production rate converges to a steady-state value R. The mapping from 
the rates to R is a concave function, so maximizing R subject to a suitable constraint on the rates is a 
convex optimization problem Il48l , IT73TI . Sensitivity analysis of the RFM with respect to the rates has been 
studied in ll49ll . These results are important in the context of optimizing the protein production rate in 
synthetic biology. Ref. lf39ll has shown that when the rates A* are time-periodic functions, with a common 
minimal period T, then every state-variable converges to a periodic solution with period T. In other 
words, the ribosomal densities entrain to periodic excitations in the rates (due e.g. to periodically-varying 
abundances of tRNA molecules). 

In ASEP with periodic boundary conditions a particle that hops from the last site returns to the first one. 
The mean field approximation of this model is called the ribosome flow model on a ring (RFMR). The 
periodic boundary conditions mean that the total number of ribosomes is conserved. Ref. 11511 analyzed 
the RFMR using the theory of monotone dynamical systems that admit a first integral. 

Both the RFM and the RFMR model mRNA translation on a single mRNA molecule. In this paper, 
we introduce a new model, called the RFM network with a pool (RFMNP), that includes a network 
of RFMs, interconnected through a dynamical pool of free ribosomes, to model and analyze simultaneous 
translation and competition for ribosomes in the cell. To the best of our knowledge, this is the first study 
of a network of RFMs. The total number of ribosomes in the RFMNP is conserved, leading to a first 
integral of the dynamics. Applying the theory of monotone dynamical systems that admit a first integral 
we prove several mathematical properties of the RFMNP: it admits a continuum of equilibrium points, 
every trajectory converges to an equilibrium point, and any two solutions emanating from initial conditions 
corresponding to an equal number of ribosomes converge to the same equilibrium point. These results 
hold for any set of rates and in particular when the RFMs in the network are not necessarily homogeneous. 
These stability results are important because they provide a rigorous framework for studying questions 
such as how does a change in one RFM affects the behavior of all the other RFMs in the network? Indeed, 
since a steady-state exists, this can be reduced to asking how does a change in one RFM in the network 
affects the steady-state behavior of the network? 

To analyze competition for ribosomes, we consider the effect of increasing one of the rates in one RFM, 
say RFM #1. This means that the ribosomes traverse RFM #1 more quickly. We show that this always 
leads to an increase in the production rate of RFM #1. All the other RFMs are always affected in the 
same manner, that is, either all the other production rates increase or they all decrease. Our analysis shows 
that this can be explained as follows. Increasing the rate A* in RFM #1 tends to increase the steady-state 
density in sites i + 1, i + 2,..., and decrease the density in site i of this RFM. The total density (i.e., the 
sum of all the densities on the different sites along RFM # 1) can either decrease or increase. In the first 
case, more ribosomes are freed to the pool, and this increases the initiation rates in all the other RFMs 
leading to higher production rates. The second case leads to the opposite result. The exact outcome of 
increasing one of the rates thus depends on the many parameter values defining the pool and the set of 
RFMs in the network. 

Our model takes into account the dynamics of the translation elongation stage, yet is still amenable 
to analysis. This allows to develop a rigorous understanding of the effect of competition for ribosomes. 
Previous studies on this topic were either based on simulations (see, for example, Ifl4l . 13) or did not 
include a dynamical model of translation elongation (see, e.g., Il23l . ||9j). For example, in an interesting 
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paper, combining mathematical modeling and biological experiments, Gyorgy et al. Il23l study the expres¬ 
sion levels of two adjacent reporter genes on a plasmid in E. coli based on measurements of fluorescence 
levels, that is, protein levels. These are of course the result of all the gene expression steps (transcription, 
translation, mRNA degradation, protein degradation) making it difficult to separately study the effect of 
competition for ribosomes or to study specifically the translation elongation step. Their analysis yields 
that the attainable output pi,p 2 of the two proteins satisfies the formula 

api + flp 2 = Y, (1) 

where Y is related to the total number of ribosomes (but also other translation factors and possibly 
additional gene expression factors), and a, fl are constants that depend on parameters such as the plasmid 
copy number, dissociation constants of the ribosomes binding to the Ribosomal Binding Site (RBS), etc. 
This equation implies that increasing the production of one protein always leads to a decrease in the 
production of the other protein (although more subtle correlations may take place via the effects on the 
constants a and (3). A similar conclusion has been derived for other models as well lfT6l Ch. 7]. 

In our model, improving the translation rate of a codon in one mRNA may either increase or decrease 
the translation rates of all other mRNAs in the cell. Indeed, the effect on the other genes depends on the 
change in the total density of ribosomes on the modified mRNA molecule, highlighting the importance 
of modeling the dynamics of the translation elongation step. We show however that when increasing the 
initiation rate in an RFM in the network, the total density in this RFM always increases and, consequently, 
the production rate in all the other RFMs decreases. This special case agrees with the results in ll23ll . 

Another recent study ||50ll showed that the hidden layer of interactions among genes arising from 
competition for shared resources can dramatically change network behavior. For example, a cascade of 
activators can behave like an effective repressor, and a repression cascade can become bistable. This agrees 
with several previous studies in the field (see, for example, lf29l . If2l). 

The remainder of this paper is organized as follows. The next section describes the new model. We 
demonstrate using several examples how it can be used to study translation at the cell level. Section [ITT] 
describes our main theoretical results, and details their biological implications. To streamline the presen¬ 
tation, the proofs of the results are placed in the Appendix. The final section concludes and describes 
several directions for further research. 

II. The model and some examples 

Since our model is based on a network of interconnected RFMs, we begin with a brief review of 
the RFM. 

A. Ribosome flow model 

The RFM models the traffic flow of ribosomes along the mRNA. The mRNA chain is divided into a 
set of n compartments or sites, where each site may correspond to a codon or a group of codons. The 
state-variable xflt), i — 1,... ,n, describes the ribosome occupancy at site i at time t, normalized such 
that xflt) = 0 [xflt) = 1] implies that site i is completely empty [full] at time t. Roughly speaking, one 
may also view xflt) as the probability that site i is occupied at time t. The dynamical equations of the 
RFM are: 


x\ = A 0 (l - xi) - Aixi(l - x 2 ), 
x 2 = Aixi(l - x 2 ) - A 2 x 2 (l - x 3 ), 
±3 = A 2 x 2 (l - x 3 ) - A 3 x 3 (l - x 4 ), 


in— 1 


in 


^n-2^n-2(l 

^71 —1*^71 — 1 (1 


Xn—l) ^n—l%n— 1 (1 
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Fig. 2. Topology of the RFM. State variable Xi(t) £ [0,1] describes the normalized ribosome occupancy level in site i at time t. The 
initiation rate is Ao, and A ; is the elongation rate between sites i and i+ 1. Production rate at time t is R(t) := A n x n {t). 


These equations describe the movement of ribosomes along the mRNA chain. The transition rates Ao,..., \ n 
are all positive numbers (units=l/time). To explain this model, consider the equation x 2 = A]Xi(l — x 2 ) — 
A 2 t 2 (1 — x 3 ). The term AiXi(l — x 2 ) represents the flow of particles from site 1 to site 2. This is 
proportional to the occupancy x\ at site 1 and also to 1 — x 2 , i.e. the flow decreases as site 2 becomes 
fuller. In particular, if x 2 = 1, i.e. site 2 is completely full, the flow from site 1 to site 2 is zero. This is 
a “soft” version of the rough exclusion principle in ASEP. Note that the maximal possible flow rate from 
site 1 to site 2 is the transition rate Ai. The term A 2 t 2 (1 — x 3 ) represents the flow of particles from site 2 
to site 3. 

The dynamical equations for the other state-variables are similar. Note that Ao controls the initiation 
rate into the chain, and that 

/?(f) . A n x n (t), 

is the rate of flow of ribosomes out of the chain, that is, the translation (or protein production) rate at 
time t. The RFM topology is depicted in Fig. [2] 

The RFM encapsulates simple exclusion and unidirectional movement along the lattice just as in ASEP. 
This is not surprising, as the RFM can be derived via a mean-field approximation of ASEP (see, e.g., J6l 
p. R345] and ll62l p. 1919]). However, the analysis of these two models is quite different as the RFM 
is a deterministic, continuous-time, synchronous model, whereas ASEP is a stochastic, discrete-type, 
asynchronous one. 

In order to study a network of interconnected RFMs, it is useful to first extend the RFM into a single¬ 
input single-output (SISO) control system: 

x\ = A 0 (l - x\)u - AiTi(1 - x 2 ), 
x 2 = AiXi(l - x 2 ) - A 2 x 2 (1 - x 3 ), 
x 3 = A 2 x 2 (1 - x 3 ) - A 3 t 3 (1 - x 4 ), 

A-1 A rl _ 2 T rl _ 2 (l X n —\) A n _iX n _i(l 

X n A n _iX n _i(l T n ) \ n x ni 

y = \ n x n . (3) 

Here the translation rate becomes the output y of the system, and the flow into site 1 is multiplied by 
a time-varying control u : M + —>- R + , representing the flow of ribosomes from the “outside world” into 
the strand which is related to the rate ribosomes diffuse to the 5’end (in eukaryotes) or the RBS (in 
prokaryotes) of the mRNA. Of course, mathematically one can absorb Ao into u, but we do not do this 
because we think of Ao as representing some local/mRNA-specific features of the mRNA sequence (e.g. 
the strength of the Kozak sequences in eukaryotes or the RBS in prokaryote). 

The set of admissible controls U is the set of bounded and measurable functions taking values in M + 
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for all time t > 0. Eq. ©, referred to as the RFM with input and output (RFMIO) |[43l , facilitates the 
study of RFMs with feedback connections. We note in passing that © is a monotone control system as 
defined in 0. From now on we write © as 

X = f(x,u), 

y = X n x n . (4) 


Let 

C n := {z G M n : z* G [0,1], i = 1, ..., n}, 

denote the closed unit cube in M n . Since the state-variables represent normalized occupancy levels, we 
always consider initial conditions x(0) G C n . It is straightforward to verify that C n is an invariant set 
of ©, i.e. for any u eU and any x(0) G C n the trajectory satisfies x(t,u ) G C n for all t > 0. 


B. RFM network with a pool 

To model competition for ribosomes in the cell, we consider a set of m > 1 RFMIOs, represent¬ 
ing m different mRNA chains. The /th RFMIO has length n, , input function u\ output function y\ and 
rates Ag, ..., X l n . The RFMIOs are interconnected through a pool of free ribosomes (i.e., ribosomes that 
are not attached to any mRNA molecule). The output of each RFMIO is fed into the pool, and the pool 
feeds the initiation locations in the mRNAs (see Fig. ©. Thus, the model includes m RFMIOs: 

x 1 — fix 1 , u 1 ), y 1 = \ l ni x ni , 


x 


= f(x m : u m ), y m = K 


X r 


and a dynamic pool of ribosomes described by 


m m 

Z = yJ ~ A o( X _ x{)Gj{z). 

j =1 3 = 1 


(5) 

( 6 ) 


where z : M + —> M describes the pool occupancy. Eq. © means that the flow into the pool is the sum 
of all output rates V J of the RFMIOs minus the total flow of ribosomes that bind to an mRNA 
molecule Kft ~ x i)Gj(z). Recall that the term (1 — x{) represents the exclusion, i.e. as the first 
site in RFMIO becomes fuller, less ribosomes can bind to it. Thus, the input to RFMIO is 


u \t) = G j( z (t))i j = 


(7) 


We assume that each Gj(-) : K + — > M + satisfies: (1) G :/ (0) = 0; Gj is continuously differentiable 
and Gj(z) > 0 for all 2 > 0 (so Gj is strictly increasing on M + ); and (3) there exists s > 0 such 
that Gj(z ) < sz for all 2 > 0 sufficiently small. These properties imply in particular that if the pool is 
empty then no ribosomes can bind to the mRNA chains, and that as the pool becomes fuller the initiation 
rates to the RFMIOs increase. 

Typical examples for functions satisfying these properties include the linear function, say, Gj(z) = z, 
and Gj(z) = aj tanh(bjz), with a 3 . b 3 > 0. In the first case, the flow of ribosomes into the first site 
of RFM is given by Aqz(1 — x\), and the product here can be justified via mass-action kinetics. The 
use of tanh may be suitable for modeling a saturating function. This is in fact a standard function in ASEP 
models with a pool ||T3l , [[ill , because it is zero when z is zero, uniformly bounded, and strictly increasing 
for z > 0. Also, for z > 0 the function tanh(z) takes values in [0,1) so it can also be interpreted as a 
probability function Il20ll . 

In the context of a shared pool, it is natural to consider the special case where G 3 (z) = G(z) for 
all j = 1,... ,m. The differences between the initiation sites in the strands are then modeled by the 
different Aq’s. 
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Fig. 3. Topology of the RFMNP. 


Note that combining the properties of G 3 with © implies that if z{0) > 0 then z(t) > 0 for all t > 0. 
Thus, the pool occupancy is always non-negative. 

Summarizing, the RFM network with a pool (RFMNP) is given by equations ©, ©, and ©. This is 
a dynamical system with d 1 + Y^Li state-variables. 

Example 1 Consider a network with m = 2 RFMIOs, the first [second] with dimension n i = 2 [ n 2 = 3]. 
Then the RFMNP is given by 

x\ = Ag(l - xDG^z) - A}^(1 - xl), 
x\ = \{x\(l-xl) - X^xl, 
xl = \l(l-xl)G 2 (z)-\lxl(l-xl), 
xl = \\x\{l-xl)-\lxl{l-xl)i 
x 3 — A 2 x 2 (1 x 3 ) — A 3 £ 3 , 

z = A ^2 + ^l x l - Ao(! - x l)Gi{z) - Aq(1 - xl)G 2 (z). (8) 

Note that this system has d = 6 state-variables. □ 

An important property of the RFMNP is, that being a closed system, the total occupancy 

m n j 

H (t) ■= z(t) + (9) 

7 = 1 i= 1 

is conserved, that is, 

H(t) = H( 0), for all t > 0. (10) 

In other words, H is a first integral of the dynamics. In particular, this means that z(t) < 11 it) = I HO) 
for all t > 0, i.e. the pool occupancy is uniformly bounded. 

The RFMNP models mRNAs that compete for ribosomes because the total number of ribosomes is 
conserved. As more ribosomes bind to the RFMs, the pool depletes, Gj(z) decreases, and the effective 
initiation rate to all the RFMs decreases. This allows to systematically address important biological 
questions on large-scale simultaneous translation under competition for ribosomes. The following examples 
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demonstrate this. We prove in Section [III] that all the state-variables in the RFMNP converge to a steady- 
state. Let e® e [0,1] denote the steady-state occupancy in site j in RFMIO #/', and let e z £ [0, oo) denote 
the steady-state occupancy in the pool. In the examples below we always consider these steady-state 
values (obtained numerically by simulating the differential equations). 

Example 2 Although we are mainly interested in modeling large-scale simultaneous translation, it is 
natural to first consider a model with a single mRNA molecule connected to a pool of ribosomes. From a 
biological perspective, this models the case where there is one gene that is highly expressed with respect 
to all other genes (e.g. an extremely highly expressed heterologous gene). 

Consider an RFMNP that includes a single RFMIO (i.e. m = 1), with dimension rii = 3, rates A] = 1, 
i = 0,1,2,3, and a pool with output function G(z ) = tanh(z). We simulated this system for the initial 
condition xj(0) = 0, z(0) = c for various values of c. Note that H(t ) = H( 0) = c. Fig. [4] depicts the 
steady-state values e 1: e 2 , e 3 of the state-variables in the RFMIO, and the steady-state pool occupancy e z . 
It may be seen that for small values of c the steady-state ribosomal densities and thus the production 
rates are very low. This is simply because there are not enough ribosomes in the network. The ribosomal 
densities increase with c. For large values of c, the output function of the pool saturates, as tanh(z) —>• 1, 
and so does the initiation rate in the RFMIO. Thus, the densities in the RFMIO saturate to the values 
corresponding to the initiation rate Ao = 1, and then all the remaining ribosomes accumulate in the pool. 
Using a different pool output function, for example G(z ) = z, leads to the same qualitative behavior, but 
with higher saturation values for the ribosomal densities in the RFMIO. (Note that the ribosomal densities 
in an RFM are finite even when A 0 —* oo m .) □ 

This simple example already demonstrates the coupling between the ribosomal pool, initiation rate, 
and elongation rates. When the ribosomal pool is small the initiation rate is low. Thus, the ribosomal 
densities on the mRNA are low and there are no interactions between ribosomes (i.e., no “traffic jams”) 
along the mRNA. The initiation rate becomes the rate limiting step of translation. On the other-hand, 
when there are many ribosomes in the pool the initiation rate increases, the elongation rates become rate 
limiting and “traffic jams” along the mRNA evolve. At some point, a further increase in the number of 
ribosomes in the pool will have a negligible effect on the production rate. 

It is known that there can be very large changes in the number of ribosomes in the cell during e.g. 
exponential growth. For example. Ref. [|8l reports changes in the range 6, 800 to 72, 000. The example 
above demonstrates how these large changes in the number of ribosomes are expected to affect the 
translational regimes; specifically, it may cause a switch between the different regimes mentioned above. 

The next example describes an RFMNP with several mRNA chains. Let l n £ M" denote the vector 
of n ones. 

Example 3 Consider an RFMNP with m = 3 RFMIOs of dimensions ni = 77.2 = n 3 = 3, and rates 

Aj = c, A? = 5, A 3 = 10, i = 0,...,3. 

In other words, every RFMIO has homogeneous rates. Suppose also that G*(z) = tanh(z), for i = 1, 2, 3. 
We simulated this RFMNP for different values of c with the initial condition z(0) = 0, £ 1 (0) = (1/2) 1 3 , 
x 2 (0) = (1/3) 1 3 and x 3 (0) = (1/4) 1 3 . Thus, H{ 0) = 3.25 in all the simulations. For each value of c, 
every state-variable in the RFMNP converges to a steady-state. Fig. [5] depicts the steady-state value e z 
and the steady-state output y l in each RFMIO. It may be seen that increasing c, i.e. increasing all the 
elongation rates in RFMIO #1 leads to an increase in the steady-state translation rates in all the RFMIOs 
in the network. Also, it leads to an increase in the steady-state occupancy of the pool. It may seem that 
this contradicts (fTOl) but this is not so. Increasing c indeed increases all the steady-state translation rates, 
but it decreases the steady-state occupancies inside each RFMIO so that the total H(t) = H{ 0) = 3.25 is 
conserved. 

Define the average steady-state occupancy (ASSO) in RFMIO by e j := 2- Ya=i e l- Fig. [6] depicts 
the ASSO in each RFMIO as a function of c. It may be seen as c increases the ASSO in RFMIO #1 
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Fig. 4. Steady-state values in the RFMNP in Example [2] as a function of the total occupancy H( 0) = c. 



Fig. 5. Steady-state outputs y l of the three RFMIOs and pool occupancy z as a function of the homogeneous transition rate c in RFMIO #1. 


decreases quickly, yet the ASSOs in the other two RFMIOs slowly increases. Since the ribosomes spend 
less time on RFMIO #1 (due to increased c) they are now available for translating the other RFMIOs, 
leading to the increased ASSO in the other mRNAs. □ 

From a biological point of view this example corresponds to a situation where accelerating one of 
the mRNA chains increases the protein production rates in all the mRNAs and also increases the number 
of free ribosomes. Surprisingly, perhaps, it also suggests that a relatively larger number of free ribosomes 
in the cell corresponds to higher protein production rates. This agrees with evolutionary, biological, and 
synthetic biology studies that have suggested that (specifically) highly expressed genes (that are transcribed 
into many mRNA molecules) undergo selection to include codons with improved elongation rates If32ll , 
If66l . If63ll . Specifically, two mechanisms by which improved codons affect translation efficiency and the 
organismal fitness are ll66l : (1) global mechanism: selection for improved codons contributes toward 
improved ribosomal recycling and global allocation; the increased number of free ribosomes improves 
the effective translation initiation rate of all genes, and thus improves global translation efficiency; and 
(2) local mechanism: the improved translation elongation rate of an mRNA contributes directly to its 
protein production rate. 
























10 



Fig. 6. Average steady-state occupancy in the three RFMIOs as a function of the homogeneous transition rate c in RFMIO #1. 


The example above demonstrates both mechanisms, as improvement of the translation elongation rates 
of one RFM increases the translation rate of this mRNA (local translation efficiency), and also of the other 
RFMs (global translation efficiency). In addition, as can be seen, the decrease in ASSO in RFMIO #1 is 
significantly higher than the increase in ASSO in the other RFMIO. Thus, the simulation also demonstrates 
that increasing the translation rate c may contribute to decreasing ribosomal collision (and possibly 
ribosomal abortion). 

We prove in Section [TIT] that when one of the rates in one of the RFMIOs increases two outcomes 
are possible: either all the production rates in the other RFMs increase (as in this example) or they all 
decrease. As discussed below, we believe that this second case is less likely to occur in endogenous genes, 
but may occur in heterologous gene expression. 

The next example describes the effect of changing the length of one RFMIO in the network. 

Example 4 Consider an RFMNP with m = 2 RFMIOs of dimensions ri\ and n 2 = 10, rates 

A- = 1 , i = 0,...,ni, 

A? = l, j = 0,...,10, 

and Gi(z) = tanh(z/200), i — 1,2. In other words, both RFMIOs have the same homogeneous rates. 
We simulated this RFMNP for different values of n\ with the initial condition z(0) = 100, x x (0) = 0 ni , 
x 2 (0) = Oio- Thus H( 0) = 100 in all the simulations. For each value of ni, every state-variable in 
the RFMNP converges to a steady-state. Fig. [7] depicts the steady-state values of z, and the steady-state 
output y l in each RFMIO. It may be seen that increasing n 1 , i.e. increasing the length of RFMIO #1 
leads to a decrease in the steady-state production rates and in the steady-state pool occupancy. This is 
reasonable, as increasing ri\ means that ribosomes that bind to the first chain remain on it for a longer 
period of time. This decreases the production rate y 1 and, by the competition for ribosomes, also decreases 
the pool occupancy and thus decreases y 2 . □ 

From a biological point of view this suggests that decreasing the length of mRNA molecules contributes 
locally and globally to improving translation efficiency. A shorter coding sequence improves the translation 
rate of the mRNA and, by competition, may also improve the translation rates in all other mRNAs. Thus, 
we should expect to see selection for shorter coding sequences, specifically in highly expressed genes and 
in organisms with large population size. Indeed, previous studies have reported that in some organisms 
the coding regions of highly expressed genes tend to be shorter lfl4l ; other studies have shown that other 
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Fig. 7. Steady-state outputs y l of the two RFMIOs and steady-state pool occupancy z/300 as a function of the length ni of RFMIO #1. 
The normalization of z by 300 is used to obtain similar magnitudes for all the values only. 



(non-coding) parts of highly expressed genes tend to be shorter ifTTl . ifTTIl . |(35l . ||64j . Decreasing the 
length of different parts of the gene should contribute to organismal fitness via improving the energetic 
cost of various gene expression steps. For example, shorter genes should improve the metabolic cost of 
synthesizing mRNA and proteins; it can also reduce the energy spent for splicing and processing of RNA 
and proteins. However, there are of course various functional and regulatory constraints that also contribute 
to shaping the gene length (see, for example OH). Our results and these previous studies suggest that 
in some cases genes are expected to undergo selection also for short coding regions, as this reduces the 
required number of translating ribosomes. 

The next section describes theoretical results on the RFMNR All the proofs are placed in the Appendix. 

III. Mathematical properties of the RFMNP 
Let 

Q . := [0, l] ni x • • • x [0, l] nm x [0, oo) 

denote the state-space of the RFMNP (recall that every x] takes values in [0,1] and z G [0, oo)). For 
an initial condition a G O, let [x(t, a) z(t, a)] denote the solution of the RFMNP at time t. It is 
straightforward to show that the solution remains in Q for all t > 0. Our first result shows that a slightly 
stronger property holds. 

A. Persistence 


Proposition 1 For any r > 0 there exists e = e{r) > 0, with e(r) —* 0 when t -g 0, such that for 
all t > t, all j — 1,..., m, all i — 1,, rtj, and all a G (O \ {0}), 

e < xl(t, a) < 1 — £, 

and 

e < z(t , a). 

In other words, after time r the solution is e-separated from the boundary of Q. This result is useful 
because on the boundary of Q, denoted f)Q, the RFMNP looses some desirable properties. For example, 
its Jacobian matrix may become reducible on dFl. Prop. [Hallows us to overcome this technical difficulty, 
as it implies that any trajectory is separated from the boundary after an arbitrarily short time. 
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B. Strong Monotonicity 

Recall that a cone K C M'" defines a partial order in M n as follows. For two vectors a, b G R n , we 
write a < b if (b — a) G K ; a < b if a < b and a ^ b\ and a -C b if (6 — a) G Int(AT). A dynamical 
system x = /(x) is called monotone if a < b implies that x(t , a) < x(t, b ) for all t > 0. In other words, 
monotonicity means that the flow preserves the partial ordering If59ll . It is called strongly monotone if a < b 
implies that x(t, a) <C x(t, b) for all t > 0. 

From here on we consider the particular case where the cone is K = M™. Then a < l> if a, < b; for 
all i, and a <C b if ai<b i for all i. A system that is monotone with respect to this partial order is called 
cooperative. 

m 

The next result analyzes the cooperativity of the RFMNR Let d := 1 + n i denote the dimension of 
the RFMNP. 


Proposition 2 For any a. b G L with a < b, 


x(t,a ) < x(t,b ) and z(t,a ) < z(t,b), 

for all t > 0. 

(ID 

Furthermore, if a < b then 



x(t, a) <C x(t, b) and z(t, a) < z(t, b), 

for all t > 0. 

(12) 


This means the following. Consider the RFMNP initiated with two initial conditions such that the 
ribosomal densities in every site and the pool corresponding to the first initial condition are smaller or 
equal to the densities in the second initial condition. Then this correspondence between the densities 
remains true for all time t > 0. 


C. Stability 
For s > 0, let 

L s := {y G Ft : l' d y = s}. 

In other words, L s is a level set of the first integral H. 

Theorem 1 Every level set L s , s > 0, contains a unique equilibrium point e/\, of the RFMNP, and for any 
initial condition a G L s , the solution of the RFMNP converges to ej Js . Furthermore, for any 0 < s < p, 

e Ls < e Lp . (13) 

In particular, this means that every trajectory converges to an an equilibrium point, representing steady- 
state ribosomal densities in the RFMIOs and the pool. Note that Proposition Q] implies that for any s > 
0, eL s G Int(Q). Eq. (fl3l) means that the continuum of equilibrium points, namely, {eL s ■ s G [0, oo)}, 
are linearly ordered. 

Example 5 Consider an RFMNP with m = 2 RFMIOs with dimensions ni = rz 2 = 1, and Gfz) = z, 
i = 1,2, i.e. 

x\ = Ag(l — x\)z — Ajxj, 
x\ = Ag(l — x\)z — X^xl, 

z = A}x} + A^x^ — Aq(1 — x\)z — Aq(1 — x\)z. (14) 

Note that even in this simple case the RFMNP is a nonlinear system. Assume that Aq = Aq = 1, and 
that A[ = Af, and denote this value simply by A. Pick an initial condition in f), and let s := x}(0) + xf(0) + 
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Fig. 8. 


Trajectories of the system in Example [5] for three initial conditions in Li. The equilibrium point GB is marked by a circle. 


z(0), so that the trajectory belong to L s for all t > 0. Any equilibrium point e = [ei e 2 e 2 ] ; e L s 
satisfies 


(1 - ei)e 2 = Aei, 
(1 - e 2 )e 2 = Ae 2 , 
ei + e 2 + e 2 = s. 


This yields two solutions 

ei = e 2 = + 2 + A — y/(s + 2 + A) 2 — 8s^ /4, (15) 

e 2 = (s - 2 - A + V(s + 2 + A) 2 -8s) /2, 

and 

ei = e 2 = + 2 + A + y/(s + 2 + A) 2 — 8s^ /4, 

e 2 = (s - 2 - A - + 2 + A) 2 - 8s) /2, 

It is straightforward to verify that in the latter solution e 2 < 0, so this is not a feasible solution. The 
solution ([f5l) does belong to L s , so the system admits a unique equilibrium in L s . Fig. [8] depicts trajectories 
of (fl4l) for three initial conditions in L\, namely, [l 0 0], [0 1 0], and [0 1/2 1/2], and the 
equilibrium point (fl5l) for s = A = 1. It may be seen that every one of these trajectories converges to e. □ 


D. Contraction 

Contraction theory is a powerful tool for analyzing nonlinear dynamical systems (see, e.g., If36l0 . with 
applications to many models from systems biology 0, lf54ll . POl . In a contractive system, the distance 
between any two trajectories decreases at an exponential rate. It is clear that the RFMNP is not a contractive 
system on 0, with respect to any norm, as it admits more than a single equilibrium point. Nevertheless, 
the next result shows that the RFMNP is non-expanding with respect to the norm: //1 = V'/ , q r \. 
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Proposition 3 For any a, b e fl, 


x(t, a) 
z(t , a) 


x(i, b ) 

z(t, b) 


< \a — 6|i, for all t > 0. 


In other words, the l\ distance between trajectories can never increase. 
Pick a6(l, and let s := l' d a. Substituting b = eL s in (fl6l) yields 


x(t, a) 
z(t, a) 


~ e L s 


<\a — e La \i, for all t > 0. 


(16) 


(17) 


This means that the convergence to the equilibrium point ex s is monotone in the sense that the (j 
distance to e^ s can never increase. Combining ([T71) with Theorem |T] implies that every equilibrium point 
of the RFMNP is semistable lf26l . 


E. Entrainment 

Many important biological processes are periodic. Examples include circadian clocks and the cell cycle 
division process. Proper functioning requires certain biological systems to follow these periodic patterns, 
i.e. to entrain to the periodic excitation. 

In the context of translation, it has been shown that both the RFM |[39l and the RFMR iBTft entrain to 
periodic translation rates, i.e. if all the transition rates are periodic time-varying functions, with a common 
(minimal) period T > 0 then each state variable converges to a periodic trajectory, with a period T. Here 
we show that the same property holds for the RFMNP. 

We say that a function / is T-periodic if f(t + T) = fit ) for all t. Assume that the Xfs in the RFMNP 
are time-varying functions satisfying: 

• there exist 0 < <5i < <5 2 such that ( t ) e [<5i, 5f\ for all t > 0 and all j e {1,, m} , i e {1,..., rij}. 

• there exists a (minimal) T > 0 such that all the Aj’s are T-periodic. 

We refer to the model in this case as the periodic ribosome flow model network with a pool (PRFMNP). 

Theorem 2 Consider the PRFMNP. Fix an arbitrary s > 0. There exists a unique function 0 S : R + —> 
Int(Q), ihat is T-periodic, and for any a £ L s the solution of the PRFMNP converges to f s . 

In other words, every level set L s of H contains a unique periodic solution, and every solution of 
the PRFMNP emanating from L s converges to this solution. Thus, the PRFMNP entrains (or phase 
locks) to the periodic excitation in the Al’s. This implies in particular that all the protein production rates 
converge to a periodic pattern with period T. 

Note that since a constant function is a periodic function for any T, Theorem [2] implies entrainment to 
a periodic trajectory in the particular case where one of the A:-\s oscillates, and all the other are constant. 
Note also that the stability result in Theorem Q] follows from Theorem [2] 

Example 6 Consider the RFMNP ([8]) with Gfz) = tanh(z), and all rates equal to one except for A|(f) = 
5 + 4sin(27rf). In other words, there is a single time-varying periodic rate in RFMIO #2. Note that all 
these rates are periodic with a common minimal period T = 1. Fig. [9] depicts the solution of this PRFMNP 
as a function of time t for 16.9 <t< 20. The initial condition is z(0) = x*(0) = 1/4 for all i,j. It may 
be seen that all the state variables converge to a periodic solution. In particular, all state variables xf(t) 
converge to a periodic solution with (minimal) period T = 1, and so does the pool occupancy z(t). 
The Xj(ty s also converge to a periodic solution, but it is not possible to tell from the figure whether there 
are small oscillations with period T = 1 or the convergence is to a constant (of course, in both cases this 
is a periodic solution with period T = 1). However, it can be shown using the first two equations in ([8]) 
that if z(t) converges to a periodic solution then so do x\(t) and x\(t). Note that the peaks in x\(t) are 
correlated with dips in x\{t), this is because when A \{t) is high on some time interval, i.e. the transition 
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Fig. 9. Trajectories of the PRFMNP in Example [6] as a function of time. 


rate from site 2 to site 3 is high, there is a high flow of ribosomes from site 2 to site 3 during this 
interval. □ 

From the biophysical point of view this means that the coupling between the mRNA molecules can 
induce periodic oscillations in all the protein production rates even when all the transition rates in the 
molecules are constant, except for a single rate in a single molecule that oscillates periodically. The 
translation rate of codons is affected among others by the tRNA supply (i.e. the intracellular abundance 
of the different tRNA species) and demand (i.e. total number of codons from each type on all the mRNA 
molecules) (see for example 1471 ). Thus, the translation rate of a codon(s) is affected by changes in the 
demand (e.g. oscillations in mRNA levels) or by changes in the supply (e.g. oscillations in tRNA levels). 
The results reported here may suggest that oscillations in the mRNA levels of some genes or in the 
concentration of some tRNA species (that occur for example during the cell cycle lf60il . |fl9l ), can induce 
oscillations in the translation rates of the rest of the genes. 

F. Competition 

We already know that any trajectory of the RFMNP converges to an equilibrium point. A natural 
question is how will a change in the parameters (that is, the transition rates) affect this equilibrium point. 
For example, if we increase some transition rate Xj in RFMIO #j, how will this affect the steady-state 
production rate in the other RFMIOs? Without loss of generality, we assume that the change is in a 
transition rate of RFMIO #1. 

Theorem 3 Consider an RFMNP with m RFMIOs with dimensions ni ,..., n m . Let A := [Aq ... A ’"' m ]' 

denote the set of all parameters of the RFMNP, and let 

e = [e! • • • ei, e? ... e 2 H2 ... ej> ... e™ ej' e (0, l)”' + " + ”“ x R ++ 

denote the equilibrium point of the RFMNP on some fixed level set of H. Pick i G {0,..., ri \}. Consider 
the RFMNP obtained by modifying Aj to A], with A! > A], Let e denote the equilibrium point in the new 
RFMNP and let e := e — e. Then 

ej < 0, and > 0, for all j G {i + 1,..., 77a}, 
sign(e*) = sign(e 2 ), for all i and all j. 

(In the case i = 0, the condition ej < 0 is vacuous.) 


(18) 

(19) 
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Increasing A] means that ribosomes flow “more easily” from site i to site / +1 in RFMIO #1. Eq. (fl8l) 
means that the effect on the density in this RFMIO is that the number of ribosomes in site i decreases, 
and the number of ribosomes in all the sites to the right of site i increases. Eq. (fl9l) describes the effect 
on the steady-state densities in all the other RFMIOs and the pool: either all these steady-state values 
increase or they all decrease. The first case agrees with the results in Example [3] above. 

Note that (fl8l) does not provide information on the change in ej, j < i. Our simulations show that any 
of these values may either increase or decrease, with the outcome depending on the various parameter 
values. Thus, the amount of information provided by (fl8l) depends on i. In particular, when A* is changed 
to A^ > A.,/ then the information provided by (fl8l) is only that 

4, < o. 

Much more information is available when i — 0. 

Corollary 1 Suppose that Aq is changed to XI > Aq. Then 

e) > 0, for all j E {1,..., rti}, (20) 

and 

e* < 0, for all i f 1 and all j , and e z < 0. (21) 

Indeed, for i — 0, (fl8l) yields ([20b . Also, we know that the changes in the densities in all other RFMIOs 
and the pool have the same sign. This sign cannot be positive, as combining this with (l20l) contradicts 
the conservation of ribosomes, so ([2Tb follows. 

In other words, increasing Aq yields an increase in all the densities in RFM #1, and a decrease in all 
the other densities. This makes sense, as increasing Aq means that it is easier for ribosomes to bind to the 
mRNA molecule. This increases the total number of ribosomes along this molecule and, by competition, 
decreases all the densities in the other molecules and the pool. Note that this special case agrees well 
with the results described in [1231 (see ([[]))• 

It is important to emphasize, however, that there are various possible intracellular mechanisms that may 
affect Xf i > 0. For example, synonymous mutation/changes (in endogenous or heterologous) genes inside 
the coding region may affect the adaptation of codons to the tRNA pool (codons that are recognized by 
tRNA with higher intracellular abundance usually tend to be translated more quickly Ifl5l ). the local folding 
of the mRNA (stronger folding tend to decrease elongation rate |[65l ), or the interaction/hybridization 
between the ribosomal RNA and the mRNA If34ll (there are nucleotides sub-sequence that tend to interact 
with the ribosomal RNA, causing transient pausing of the ribosome, and delay the translation elongation 
rate). Non synonymous mutation/changes inside the coding region may also affect the elongation for 
example via the interaction between the nascent peptide and the exit tunnel of the ribosome ff37l . Il55ll . 
In addition, intracellular changes in various translation factors (e.g. tRNA levels, translation elongation 
factors, concentrations of amino acids, concentrations of Aminoacyl tRNA synthetase) and, as explained 
above, the mRNA levels can also affect elongation rates. Furthermore, various recent studies have demon¬ 
strated that manipulating the codons of a heterologous gene tend to result in significant changes in the 
translation rates and protein levels of the gene If32ll . |[72]| , [|5]| . 

Thus, our study is relevant to fundamental biological phenomena that are not covered in models that 
do not take into account the elongation dynamics. 

Example 7 Consider the RFMNP in ([8b with Gfz ) = z, Aj = Aq = 1, A^ = A 2 = 0.1, A| = 1, and initial 
condition (1/4) 1 6 . We consider a range of values for A3,. For each fixed value, we simulated the dynamics 
until steady-state for two cases: A} = 1 and A} = 10. Fig. [TOl depicts e\,e\ for the various fixed values 
of A3,. It may be seen that we always have e\ < 0 and e\ > 0. Fig. QT] depicts ef i = 1, 2, 3, and e z for 
the various fixed values of A3,. It may be seen that for a small value of A 2 all the ef’s and e z are negative, 
whereas for large values of A3, they all become positive. □ 
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Fig. 10. Steady-state perturbations in RFM in the RFMNP in Example [7] when changing A) = 1 to A) = 10 for various values of A 2 . 



Fig. 11. Steady-state perturbations in RFM #2 in the RFMNP in Example [7] when changing A) = 1 to A) = 10 for various values of 


Theorem [3] implies that when the codons of a gene are modified into “faster codons” (either via synthetic 
engineering or during evolution) then either all the translation rates of the other genes increase, or they all 
decrease. However, Theorem [3] does not provide information on when each of these two cases happens. 
In order to address this, we need to calculate derivatives of the equilibrium point coordinates with respect 
to the rates. The next result shows that these derivatives are well-defined. Denote the mapping from the 
parameters to the unique equilibrium point in Int(fi) by a, that is, e* = a*-(A, H(0)), i = 1 
j = l,...,ni. 

Proposition 4 The derivative (A. 11(0) ) exists for all i,j,p,q. 

The next example uses these derivatives to obtain information on the two cases that can take place as 
we change one of the rates. 

Example 8 Consider an RFMNP with m = 2 RFMIO’s with lengths n and l. To simplify the notation, 
let e = [ei,..., e n ] \v = [fi,..., e^] ] denote the equilibrium point of RFMIO #1 [RFMIO #2], and 
let Aj, i — 0,..., n, denote the rates along RFMIO #1. Suppose that Ai is changed to Ai. Differentiating 
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the steady-state equations 

AoGi(e z )(l ei) Aiei(l e 2 ) ••• A n _ie n _i(l c n ) A n e n , 

n i 

^ v j + e *=-^(o) ; 

*=i j=i 

w.r.t. Ai yields 

^oG^e^e^l — e i) — ^oG < i(e,)e , 1 

e 'z + J2 v 'j 

3 = 1 

where we use the notation f := -£^f. These two equations yield 

l n —1 

(A n + AoG f, 1 (e^)(l — ei))e^, + A n = (AoGi(e z ) — A n )e , 1 — X n e^. 

j=i *=2 

Recall that Gi(e z ) > 0, G\ (e z ) > 0, A j > 0 for all j, and 0 < e p < 1 for all p. Also, by Theorem [ 3 ] 

e\ < 0, e'j > 0, for all j e 2, and sign(e^) = sign(i^), for all k. Thus, 

sign«) = sign(u') = sign f (AoG^e*) - A n )e' 1 - A n e'J . (22) 

This means that the sign of the change in the densities in all the other RFMIO’s and the pool depends 

on several steady-state quantities including terms related to the initiation rate A 0 G| (e 2 ) and exit rate A n 
in RFMIO #1, and also the change in the total density Y^iZi e \ i n this RFMIO. 

In the particular case n — 2 (i.e., a very short RFM), Eq. (l22l) becomes: 

sign(e^) = sign(n') = sign(A 2 - A 0 Gi(e z )). (23) 

Note that A 0 Gi(e 2 ) [A 2 ] is the steady-state initiation [exit] rate in RFMIO #1. Thus, A 2 — A 0 Gi (e z ) > 0 
means that it is “easier” for ribosomes to exit than to enter RFMIO #1, and in this case (1231) means that 
when Ai is increased the change in all other densities will be positive. This is intuitive, as more ribosomes 
will exit the modified molecule and this will improve the production rates in the other molecules. On the 
other-hand, if A 2 — A 0 Gi (ez) < 0 then it is “easier” for ribosomes to enter than to exit RFMIO #1, so 
increasing Ai will lead to an increased number of ribosomes in RFMIO #1 and, by competition, to a 
decrease in the production rate in all the other RFMIOs. □ 

Note that in the example above increasing Ai always increases the steady-state production rate R = A 2 e 2 
in RFM #1 (recall that e' 2 > 0). One may expect that this will always lead to an increase in the production 
rate in the second RFMIO as well. However, the behavior in the RFMNP is more complicated because 
the shared pool generates a feedback connection between the RFMIO’s in the network. In particular, the 
effect on the other RFMIO’s depends not only on the modified production rate of RFMIO #1, but on 
other factors including the change in the total ribosome density in RFMIO #1 (see (l22l) f. 

This analysis of a very short RFM suggests that the steady-state initiation rate of the mRNA with the 
modified codon plays an important role in determining the effect of modifications in the network. If this 
initiation rate is relatively low (so it becomes the rate limiting factor), as believed to be the case in most 
endogenous genes ff28l . then the increase in the rate of one codon of the mRNA increases the translation 
rate in all the other mRNAs, whereas when this initiation rate is high then the opposite effect is obtained. 
This latter case may occur for example when a heterologous gene is highly expressed and thus “consumes” 
some of the available elongation/termination factors making the elongation rates the rate limiting factors. 



i— 1 
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IV. Discussion 

We introduced a new model, the RFMNP, for large-scale simultaneous translation and competition for 
ribosomes that combines several RFMIOs interconnected via a dynamic pool of free ribosomes. To the best 
of our knowledge, this is the first model of a network composed of interconnected RFMIOs. The RFMNP 
is amenable to analysis because it is a monotone dynamical system that admits a non-trivial first integral. 
Our results show that the RFMNP has several nice properties: it is an irreducible cooperative dynamical 
system admitting a continuum of linearly ordered equilibrium points, and every trajectory converges to an 
equilibrium point. The RFMNP is also on the “verge of contraction” with respect to the l\ norm, and it 
entrains to periodic transition rates with a common period. The fact that the total number of ribosomes in 
the network is conserved means that local properties of any mRNA molecule (e.g., the abundance of the 
corresponding tRNA molecules) affects its own translation rate, and via competition, also globally affects 
the translation rates of all the other mRNAs in the network. 

An important implication of our analysis and simulation results is that there are regimes and parameter 
values where there is a strong coupling between the different “translation components” (ribosomes and 
mRNAs) in the cell. Such regimes cannot be studied using models for translation of a single isolated 
mRNA molecule. The RFMNP is specifically important when studying highly expressed genes with many 
mRNA molecules and ribosomes translating them because the dynamics of such genes strongly affects 
the ribosomal pool. For example, changes in the translation dynamics of a heterologous gene which is 
expressed with a very strong promoter, resulting in very high mRNA copy number should affect the 
entire tRNA pool, and thus the translation of other endogenous genes. Highly expressed endogenous 
genes ’’consume” many ribosomes. Thus, a mutation that affects their (“local”) translation rate is expected 
to affect also the translation dynamics of other mRNA molecules. Studying the evolution of such genes 
should be based on understanding the global effect of such mutations using a computational model such 
as the RFMNP. 

On the other hand, we can approximate the dynamics of genes that are not highly expressed (e.g., a 
gene with mRNA levels that are 0.01% of the mRNA levels in the cell) using a single RFM. In this case, 

the relative effect of the mRNA on all other mRNAs is expected to be limited. 

Our analysis shows that increasing the translation initiation rate of a heterologous gene will always 
have a negative effect on the translation rate of other genes (i.e. their translation rates decrease) and vice 
versa. The effect of increasing [decreasing] the translation rate of a codon of the heterologous gene on the 
translation rate of other genes is more complicated: while it always increases [decreases] the translation 
rate of the heterologous gene it may either increase or decrease the translation rate of all other genes. The 
specific outcome of such a manipulation can be studied using the RFMNP with the parameters based on 
the heterologous genes and the host genome. 

Our analysis suggests that the effect of improving the transition rate of a codon in an mRNA molecule 
on the production rate of other genes and the pool of ribosomes depends on the initiation rate in the 

modified mRNA. When the initiation rate is very low the effect is expected to be positive (all other 

production rates increase). However, if the initiation rate is high the effect may be negative. This may 
partially explain the selection for slower codons in highly expressed genes that practically decrease the 
initiation rate lf67l . If63ll . This relation may also suggest a new factor that contributes to the evolution 
of highly expressed genes towards higher elongation and termination rates (i.e., the tendency of highly 
expressed genes to include “fast” codons). Indeed, lower elongation rates (and thus a relatively high 
initiation rate) may decrease the production rates of other mRNAs that are needed for proper functioning 
of the organism. 

The RFM, and thus also the model described here, do not capture certain aspects of mRNA translation. 
For example, eukaryotic ribosomes may translate mRNAs in multiple cycles before entering the free ribo¬ 
somal pool IfTOll . [[43ll . OTI . This phenomenon may perhaps be modeled by adding positive feedback 1431 
in the RFMNR In addition, different genes are transcribed at different rates, resulting in a different number 
of (identical) mRNA copies for different genes. This can be modeled using a set of identical RFMs for 
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each gene. Such a model can help in understanding how changes in mRNA levels of one gene affect the 
translation rates of all the mRNAs. The analysis here suggests that modifying the mRNA levels of a gene 
will affect the translation rates of all other genes in the same way. These and other aspects of biological 
translation may be integrated in our model in future studies. Ref. Il25ll develops the notion of the realizable 
region for steady-state gene expression under resource limitations, and methods for mitigating the effects 
of ribosome competition. Another interesting research direction is studying these topics in the context of 
the RFMNP. 

The results reported here can be studied experimentally by designing and expressing a library of 
heterologous genes f72ll . If5l , |[32l . The effect of the manipulation of a codon (i.e, increasing or decreasing 
its rate) of the heterologous gene on the ribosomal densities and translation rates of all the mRNAs 
(endogenous and heterologous) can be performed via ribosome profiling |[27il in addition to measurements 
of mRNA levels, translation rates, and protein levels lf56fl . 

We believe that networks of interconnected RFMIOs may prove to be powerful modeling and analysis 
tools for other natural and artificial systems as well. These include communication networks, intracellular 
trafficking in the cell, coordination of large groups of organisms (e.g., ants), traffic control, and more. 
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Appendix: Proofs 

Proof of Proposition]]} We require the following result on repelling boundaries and persistence. 
Lemma 1 Consider a time-varying system 

x = f(t,x ) (24) 

whose trajectories evolve on Q := I\ x. f 2 x. ... x I n C M", where each f is an interval of the form 
[0, A], A > 0, or [0, oo). Suppose that the time-dependent vector field f = [/i,..., / n ] has the following 
two properties: 

• the cyclic boundary-repelling property (CBR). For each 5 > 0 and each sufficiently small A > 0, 
there exists K = K (<5, A) > 0 such that, for each k — 1, n and each t > 0, the condition 

Xkif) A A, and Xk-iif) > S (25) 

(all indexes are modulo n) implies that 

f k (t,x)>K ; (26) 

• for any i G (1,..., n}, and any s > 0, xfis) > 0 implies that 

xfit) > 0, for all t > s. (27) 

Then given any r > 0 there exists £ = s(t) > 0, with e{r) —> 0 as r —* 0, such that every solution x(t, a), 

with a f 0, satisfies 

xfit, a) > e for all i G {1,..., n} and all t > r. 

In other words, the conclusion is that after an arbitrarily short time every xft, a) is separated away from 
zero. 

Proof of Lemma Q] Pick r > 0 and a f 0. Then there exists i E {1,... ,n} such that a* > 0. Since the 
(CBR) condition is cyclic, we may assume without loss of generality that x n (0) > 0. Then (l27l) implies 
that there exists <5 > 0 such that x n (t) > <5 for all t E [0, r/n]. 

Fix A > 0 such that (CBR) holds. Let K = K(S,A), and define £\ := min{A, Kr/n}. Let t 0 6 

[0,r/n] be such that x\(t 0 ) > ey. Such a t 0 exists, since by property (CBR), x\(t) < ey < A for all 
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t E [0 ,r/n] would imply that x 1 (t) = fi(t,x(t)) > K > 0 for all t E [0, r/n], which in turn implies 
Xi(r/n) > Xi(0) + Kr/n > Kr/n > e i, contradicting xi(r/n) < ey. We claim that also xi(i) > sy for 
every t >t 0 . Indeed, suppose otherwise. Then, there is some t\ > to such that £ := xi(fi) < ey. Let 

cr := minjf > t 0 \ ay (t) < £} > t 0 . 

As ay(cr) < £ < ey < A, and x n (a) > 0 property (CBR) says that iy(cr) = fi(a,x(a)) > 0, so it follows 
that x\ (t) > 0 on an interval [a — c, cr], for some c > 0. But then ay (cr — c) < ay (a) < £, which contradicts 
the minimality of a. Thus X\ (t) > £\ for all t > t 0 , and in particular 

x\(t) > £i, for all t > r/n. (28) 


Let K ■= K(£i, A), and define £ 2 := min{A, Kr/n}. Let t 0 E [r/n, 2r/n\ be such that x 2 (fo) A £ 2 • 
Such a to exists, since by property (CBR) and (l28l) . x 2 (t) < £2 < A for all t E [r/n, 2r/n\ would 
imply that x 2 (t) = f 2 (t,x(t)) > K > 0 for all t E [r/n,2r/n], which in turn implies x 2 (2 r/n) > 
X 2 (r/n) + Kr/n > Kr/n > £ 2 , contradicting x 2 (2r/n) < £ 2 . We claim that also x 2 (t) > £ 2 for every 
t > t 0 . Indeed, suppose otherwise. Then, there is some t x > t 0 such that £ := x 2 (ti) < e 2 - Let 

cr := minjf > t 0 \ x 2 (t) < £} > fo- 

As x 2 (a) < £ < £ 2 < A, and xi(cr) > £\ property (CBR) says that x 2 (a) = f 2 (a,x(a)) > 0, so it 
follows that x 2 (t) > 0 on an interval [cr — c, cr], for some c > 0. But then x 2 (a — c) < x 2 (a) < £, which 
contradicts the minimality of a. Thus x 2 (t) > £ 2 for all t > to, and in particular 

x 2 (t) > e 2 , for all t > 2 r/n. 

Continuing in this manner we have that for every i E {1,..., n} there exists c, > 0 such that 

Xi(t) > £i, for all t > irIn. 


Thus, 


Xi(t) > min{£!,..., e n }, for all t > r and alii = 1,..., n, 

and this completes the proof of Lemma Q] □ 

We can now prove Proposition Q] Consider the case m — 1, i.e. the RFMNP is given by 

xi = A 0 (l - xi)G(x n+ i) - Aixi(l - x 2 ), 
x 2 = AiXi(l - x 2 ) - \ 2 x 2 (l - x 3 ), 


X n 1 A n _ 2 x n _ 2 (l X n —\) \ n —iX n —i(\ Xn), 

x n A n _ix n _i(l x n ) X n x n , 

Xu • 1 X n x n A 0 (l x\)G(x n -\-i) , (29) 

where we write x n+ \ instead of z. The proof in the case where m > 1 is similar. We begin by showing 
that (l29l) satisfies the properties in Lemma Q] on = [0,1]" x [0, 00 ). Fix an arbitrary 5 > 0. If x\ < A 
and x n+ i > 5 then 

fi = A 0 (l - xi)G(x n+ i) - AiXi(l - x 2 ) 

>K U 

where K x ■= A 0 (l — A )G(8) — AiA. Now pick 1 < k < n. If Xfc < A and > 5 then 

fk Xfc) AfcX^(l Xfc_|_x) 

> Kk, 
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where 

Kk '■= Afc_i5(l — A) — A fc A. 

If x n < A and x n _i ><5forl<i<n — 1 then 

fn A n _iX n _i(l X n ^) \ n X n 

> K n -\. 

Finally, if x n+ 1 < A and x n > S then 

fn+l X n X n A 0 (l X\)G(x n - |_i) 

where K n+ 1 := A n 5 — A 0 G(A). Thus, (l26l) holds for K := minjA'i,..., K n+ ±}, and clearly K > 0 for 
all A > 0 sufficiency small. Thus, ([29b satisfies (CBR). 


To show that (l29l) also satisfies (f27l) . note that for all k G {1,..., n} and all x e fl, x k > —X k x k , and 
that by the properties of G(-), x n+ i > —sx n+ 1 for all x n+ i > 0 sufficiently small. Thus, (l29l) satisfies all 
the conditions in Lemma [Hand this implies that for any r > 0 there exists e = e{r j 2) > 0 such that 

Xi(t, a) > e, for all? e {1,..., n + 1} and all t > t/2. (30) 

This proves the first part of Proposition HI To complete the proof, define 

Ui(t) := 1 - x n+ i_i(t), i = (31) 

and y n +\(t) := x n+1 (t). Then 

Vi = A n (l - Vi) - A n _ij/i(l - y 2 ), 
y 2 = A n _iyi (1 — y 2 ) — A n _ 2 |/2 (1 — 2/3); 

i)n— 1 X 2 y n — 2(1 Vn—l) Aiy n _i(l yri)i 

ijn Aiy n _i(l y n ) AoG(y n _|_i)y n , 

ijn+l A n (l yi) A 0 G(yn+l)yn. 

The first n equations here are an RFM with a time-varying exit rate A 0 G(y, H -i (t)). We already know 
that y n+ i(f) > e for all t > r/2 and the results in lf39ll imply that there exists e 1 = E\ (r) > 0 such that 

yi{t, a) > £ 1 , for all* e {1,, n} and all t > r. 

Combining this with OTI) completes the proof of Prop. Q] □ 


Proof of Proposition^ The Jacobian matrix J of the RFMNP has the form 


Jl,l 

0 

0 


0 

J 2 , 2 

0 

^2,m+l 

0 

0 

Jm,m 

J m,m+l 


Jm+ 1,2 

J m+l,ra 

J ra+l,m+l_ 


( 32 ) 
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Here J i}i G M niXn \ i = 1,..., m, is the Jacobian of RFMIO jfi given by 


xiGi(z) - \\(i ~ xi) 
Al(i-4) 


Ajxl 0 

- X 2 (l-x\) ^2 X 2 

y 2 (l ~ X V) ^ 2 x 2 — A|(l — x\) 


0 

0 

0 


0 

0 

0 


0 


0 

0 


0 

0 


0 

0 



and the other blocks are 



i — 1,..., m, and 


m 


Jm+l,m+l ~ ~ Ajj(l — x{)G'j(z) G M. 


Since x* G [0,1] and A* > 0, every off-diagonal entry of ./,,, is non-negative. Since z > 0, every 
entry of J m +i,j, Ji,m+u i = 1, is also nonnegative. We conclude that every non-diagonal entry 

of J is non-negative, and this implies (fill) (see 11591 1. To prove (fl2l) . note that for any point in Int(Q) 
(i.e., Xj G (0,1), j — 1,..., d — 1, and z > 0) every entry on the super- and sub-diagonal of J iti is strictly 
positive. Also, Gi(z) > 0, i — 1,..., m. This implies that the matrix J in (l32l) is irreducible. Combining 
this with Prop. Q] completes the proof. □ 

Proof of Theorem [0 Since the RFMNP is a cooperative irreducible system on Int(Q) with a non-trivial 
first integral, Thm. Q] follows from combining Prop. Q] with the results in [j46l (see also ||5T1 . P31 and lf33l 
for related ideas). □ 

Proof of Proposition^ J Recall that the matrix measure //1 (•) : M nxn —> M induced by the norm is 


pi(A) = max{ci(A), ...,c n (A)}, 


where q(A) := an + ^2^ \aki\, i.e. the sum of entries in column i of A, with the off-diagonal entries 
taken with absolute value |[68ll . For the Jacobian of the RFMNP (l32l) . we have c l (-f(x)) = 0 for all i 
and all i G fl, so pi(J(x)) = 0. Now (fT6l) follows from standard results in contraction theory (see, 

e.g., El). □ 

Proof of Theorem C3 Write the PRFMNP as x = f(t,x). Then f(t,y) = fit + T, y) for all t and y. 
Furthermore, AT in ([9]) is a non trivial first integral of the dynamics. Now Theorem [2] follows from the 
results in [l6Tll (see also l[30ll ). The fact that (fr s G Int(O) follows from Proposition Q] □ 

Proof of Theorem Q1 To simplify the presentation, we prove the theorem for the case m — 2 and change 
some of the notation. The proof in the general case is very similar. Let n [(] denote the dimension of 
the first [second] RFMIO, and let A* [(,] denote the rates in the first [second] RFMIO. We denote the 
state-variables of RFMIO # 1 by x % , i = 1 ,,n, and those of the second RFMIO by y t , i = 1, 

Let en i = 1 ,... ,n, [vj, j = 1,..., f\ denote the equilibrium point of the first [second] RFMIO. Then 
we need to prove that 


e-i < 0, and e 3 > 0, for all j G {i + 1,..., n}, 
sign(Fi) = • • • = sign(^) = sign(e~), 


(33) 

(34) 
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where iij Vj — v r At steady state, the RFMNP equations yield: 

AoGi(e 2 )(l — ei) = Aiei(l — e 2 ), 

= A 2 e 2 (l — e 3 ), 

= Ase 3 (l — e 4 ), 

An—iTn—i(f e n ) 

An^ni (35) 

CoG 2 (e 2 )(l - vi) = Ciui(l - v 2 ), 

= C 2 ^ 2(1 - V 3 ), 

= ( 3 ^ 3(1 - v 4 ), 

= 0-i^-i(i - 

= 0^, (36) 

A n G n + = AoG'i(e 2 )(l — ei) + C 06 * 2 ( 63 )(1 — vi). (37) 

Also, since H is a first integral 

n i 

£ e *+E Vj + e z = H( 0). (38) 

*=1 j=i 

Pick i 6 {1,..., n — 1}. Consider a new RFMNP obtained by modifying A, in RFMIO #1 to A*, 
with Aj > A, : . Then the equations for the modified equilibrium point are: 

AoGi(e^)(l — e 3 ) = Aiei(l — e 2 ), 

= A 2 e 2 (l — e 3 ), 

AjCj(l 6j+l)j 

A n _ie n _i(l e n ) 

A rfini 

CoG , 2 (e z )(l - vi) = Ciui(l - v 2 ), 

= ( 2 ^ 2(1 - U 3 ), 

= G-iG-i(l - vi) 

= CtVi, 

A nV-n + CeVi = \oGi{e z ){l — e 3 ) + (oG 2 (e z )(l — vi) 

Since the initial condition remains the same, 

n i n i 

ei + vj +e z = ^^^= 77 ( 0 ), 

Z— 1 J = 1 Z=1 J = 1 


(39) 


(40) 

(41) 
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SO 


n l 

'^y ^ Cj “I - y 'y ^ Uj ~\~ & z o. 

*=i j=i 


( 42 ) 


The last equality in (1361) yields 


Q-i v e-\ — 


Cm 

1 - Vi 


The right-hand side here is increasing in vg, and the left-hand side is increasing in vg-i (recall that Qt 
and Ct-i are the same for both the original and the modified RFMNP), so a change in \ must lead 
to sign(u£_ 1 ) = sign(fy). Using (1361) again yields 


> Cm 

U-2V(-2 — - -, 

1 - Vg-! 


so sign({y_ 2 ) = sign(£y_i) = sign(f^). Continuing in this way yields 


signal) = sign(n 2 ) = • • • = sign ({>*). 


(43) 


By d36]), 


G 


2 eri = 


Cm 


Co(l — Wl) " 


Since Gi(p) is strictly increasing in p, combining this with (1431) implies that sign(e 2 ) = sign(i) 1 ) = 
sign(n 2 ) = • • • = sign(t^). This proves (1341) . To prove (l33l) . note that arguing as above using (l39l) yields 


sign(e„) = sign(e n _i) = • • • = sign(e i+ i). 

Seeking a contradiction, assume that sign(ej) = sign(e n ). By (l39l) . 

_ ^ii&n 

Ui^) 

so sign(ej-i) = sign(e n ), and continuing in this fashion yields 


sign(e n ) = sign(e n _i) = • • • = sign(ei). 


(44) 


By (EH), 


Gi(e z ) — 




Aq( 1 — 6l) 


and combining this with (1441) implies that sign(ei) = ... = sign(e„) = sign(e~). Since also sign(u 1 ) = ... = 
sign(t^) = sign(e 2 ), it follows that all the differences have the same sign, and this contradicts (1421) . We 
conclude that if e* ^ 0 then sign(gj) ^ sign(e i+l ) = sign(e i+2 ) = • • • = sign(e„). To complete the proof 
of (1331) . we need to show that e, < 0. Seeking a contradiction, assume that e* > 0, so e i+ i < 0,..., e n < 0. 
Thus, 

Cj+1 — ('i- 1-1) ■ ■ ■ ■ ^ e n . 


By (1351) . A* 


e »(l— e x+l) ’ 


SO 


1 i A i e n e n 

Ar?, ^i(l Cj+l) 6j( 1 Cj-l-l) 

< 0. 


This contradiction completes the proof for the case where i e {1,... ,n — 1}. The proof for the case i — 0 
and i = n is similar. □ 
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Proof of Proposition 0 In Il45ll . Mierczynski considered an irreducible cooperative dynamical sys¬ 
tem, x = f(x), that admits a non trivial first integral II(x) with a positive gradient. Let S x := {v G 
R n : V II r (x)v = 0}, and consider the extended system 

x = f(x), 

Sx = J(x)Sx, 

where J := is the Jacobian, with initial condition x(0) = x 0 , Sx( 0) = Sx 0 G S x \ {0} . Mierczynski 
shows that there exists a norm, that depends on x, such that 

\Sx(t)\ x ( t ) < l&coL, for a111 > °- 

(For a general treatment on using Lyapunov-Finsler functions in contraction theory, see At the 

unique equilibrium point e this yields 

I exp(J(e)t)<Lc 0 |e < |&Eo|ej f° r a U t > 0- 

This implies that the matrix obtained by restricting J(e) to the integral manifold is Hurwitz and thus 
nonsingular. Invoking the implicit function theorem implies that the mapping from A to e can be identified 
with a C 1 function. □ 
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